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Abstract. Due to their algorithmic simplicity and high accuracy, force-based model coupling 
techniques are popular tools in computational physics. For example, the force-based quasicontinuum 
approximation is the only known pointwise consistent quasicontinuum (QC) approximation for 
coupling a general atomistic model with a finite element continuum model. In this paper, we present 
a detailed stability and error analysis of this method. Our optimal order error estimates provide 
a theoretical justification for the high accuracy of the force-based QC approximation: they clearly 
demonstrate that the computational efficiency of continuum modeling can be utilized without a 
significant loss of accuracy if defects are captured in the atomistic region. 

The main challenge we need to overcome is the fact (which we prove) that the linearized QC 
operator is typically not positive definite. Moreover, we prove that no uniform inf-sup stability 
condition holds for discrete versions of the W 1 ' p -W 1 ' 9 "duality pairing" with 1/p + 1/q = 1, if 
1 < p < oo. We therefore derive an inf-sup stability condition for a discrete version of the W 
W "duality pairing" which then leads to optimal order error estimates in a discrete W 1 



1. Introduction 

Localized defects in materials typically interact with elastic fields far beyond the defects' atomic 
neighborhood. Accurately computing the structure of localized defects requires the use of atomistic 
models; however, atomistic models are too computationally demanding to be utilized for the entire 
interacting system. The goal of atomistic to continuum coupling methods such as the quasicon- 
tinuum (QC) method is to use the computationally intensive, fully atomistic calculations only in 
regions with highly non-uniform deformations such as neighborhoods of dislocations, crack tips, and 
grain boundaries; and to use a (local) continuum model in regions with nearly uniform deformations 
to reduce the number of degrees of freedom. 

The initial computational results obtained with the QC method have excited the materials sci- 
ence community with the promise of the simulation of heretofore inaccessible multiscale materials 
problems [19,26,33]. Variants of the QC method have continued to be developed with the intro- 
duction of adaptive methods, improved mesh generation, and faster solvers [1,2,8,15,24,29,30]; 
yet, in common with many other multiscale methods it lacks the theoretical basis needed to give 
predictive computational results. 

During the past few years, a mathematical structure has been given to the description and anal- 
ysis of various flavors of the QC method, clarifying the relation between different approximations 
and the corresponding sources of error [3, 7, 10, 12-14, 16, 17, 22, 23, 25, 28]. In the present paper, 



Date: March 8, 2009. 

2000 Mathematics Subject Classification. 65Z05,70C20. 

Key words and phrases, quasicontinuum, force-based, atomistic to continuum, stability, coercivity, error estimates. 
This work was supported in part by DMS-0757355, DMS-0811039, the Department of Energy under Award 
Number DE-FG02-05ER25706, the Institute for Mathematics and Its Applications, the University of Minnesota 
Supcrcomputing Institute, the University of Minnesota Doctoral Dissertation Fellowship, and the EPSRC critical 
mass programme "New Frontier in the Mathematics of Solids." 



STABILITY, INSTABILITY, AND ERROR OF THE FORCE-BASED QC APPROXIMATION 



2 



we contribute to this effort by providing a detailed stability and error analysis of the force-based 
quasicontinuum approximation. 

Considerable concern has been generated by the discovery that early QC approximations exhibit 
"ghost" forces in the atomistic to continuum interface when the material is subject to a uniform 
strain, that is, they do not satisfy the "patch test" criterion of computational mechanics. The first 
remedy, which is still commonly employed, is known as the ghost force correction. It first applies 
a dead load that corrects the ghost forces at the current state of a continuation process, then 
increments the parameter value for the process, then reminimizes the energy at the new parameter 
value and dead load, and finally recomputes the ghost force corrections for use as a dead load at 
the next step of the continuation process [26,30]. In [7,8] this process was identified as an iterative 
method to approximate the solution to the equilibrium equations for a purely force-based coupling 
approach, which we label the force-based quasicontinuum ( QCF) approximation. This formulation of 
the QCF method has enabled the development of more efficient iterative and continuation methods 
for its solution and a more precise understanding of the error [7, 8] . Related force-based modeling 
approaches, which couple an atomistic region with a continuum region modeled by linear elasticity 
can be found in [20,31]. 

Other research groups have proposed QC approximations that utilize special interfacial atoms 
at the atomistic to continuum interface in an attempt to develop an energy-based QC method that 
does not suffer from the ghost force problem mentioned above [12,32]. The quasi-nonlocal (QNL) 
approach [32] is easy to implement and removes ghost forces for short range interactions (depending 
on the lattice structure), but ghost forces remain for longer range interactions. This method was 
generalized in the reconstruction approach [12] which, in theory, allows for the elimination of 
all ghost forces; however, explicit methods have only been constructed for planar interfaces so 
far. Moreover, a computationally efficient implementation of this method, that can be used with 
adaptive atomistic to continuum modeling algorithms, has yet to be proposed. 

Both of the above methods [12, 32] couple the original atomistic model to a new atomistic 
model with local interactions. To allow for the reduction of degrees of freedom by piecewise linear 
interpolation in the continuum region as in the finite element method, it is necessary to further 
couple this local atomistic model to a volume-based local model. However, it is not known how to 
couple a local atomistic model to a volume-based local model along a nonplanar interface without 
introducing ghost forces [12]. In contrast, the force-based quasicontinuum approximation allows 
arbitrary atomistic to continuum interfaces and coarsening without ghost forces. 

Rather than computing forces from a total energy, the force-based quasicontinuum approximation 
directly assigns forces using a simple rule: the force on an atom in the atomistic region is computed 
from the force law of the atomistic model, while the force on a degree of freedom in the continuum 
region is computed from the force law of the continuum (finite element) approximation [7,8]. There 
is no modification of these equations near the atomistic to continuum interface and it is therefore 
easy to see that the QCF equilibrium equations are satisfied exactly by a material under uniform 
strain, that is, there are no ghost forces in this approximation. Moreover, we will show below that 
the QCF approximation has an 0(e 2 ) truncation error in the atomistic to continuum interface for 
all smoothly varying strains. By contrast, it has been shown in [9] that, even when it succeeds in 
removing ghost forces, the QNL method has an 0(1) truncation error in the atomistic to continuum 
interfaces for a nonuniform but smooth strain. (This is nevertheless a significant improvement over 
the 0(l/e) truncation error in the original QC method.) 

With the exception of [28], error analyses of energy-based QC methods have utilized the coercivity 
(positive-definiteness) of the linearization of the quasicontinuum equilibrium equations about the 
energy-minimizing solution [9,13,14]. A recent attempt to establish an error analysis for the QCF 
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method has presented an invalid proof of coercivity of the linearized equilibrium equations and an 
error analysis based on this incorrect coercivity result [27]. In the present paper we prove that, 
typically, the linearization of the QCF equilibrium equations is not coercive (cf. Theorem 14. lj) . and 
consequently, our error analysis will be based on a more general inf-sup stability condition. However, 
even this more general approach will fail unless one chooses the norms particularly carefully. We 
show in Section [5] that the linearized QCF operator is stable with respect to a discrete version of 
the W*->°°-W ' pairing, uniformly with respect to the number of atoms, but we show in Section [7] 
that it is not uniformly stable with respect to any other W 1,p -W 1,q pairing where - + ~ = 1. 

Our goal in this paper is to clearly present our techniques in the simplest setting. For this reason, 
we restrict our presentation to a one-dimensional chain of atoms which interact with nearest and 
next-nearest neighbors. To further simplify the setting, we consider a linearization of the force-based 
equilibrium equations about a uniform strain. Although the QCF approximation can be directly 
formulated and implemented with mesh coarsening in the continuum region, we only consider the 
modeling error due to the QCF approximation itself and do not consider the coarsening error. 
Each of these extensions deserve a careful analysis in order to firmly establish the mathematical 
foundation of the force-based quasicontinuum approximation. 

The main result of the present paper is that the strain error for the QCF approximation is 0(e 2 ), 
where e is the lattice spacing scaled by the material dimension. The prefactor for the e 2 error term 
is a maximum norm of a third divided difference of the displacement restricted to the continuum 
region only. Thus, our analysis predicts the observed high accuracy of the QCF method when 
defects are modeled in the atomistic region. 

In Section [21 we present a detailed description of the force-based quasicontinuum approximation 
and a first estimate of the truncation error with respect to the fully atomistic model. In Section El 
we show how to formulate the QCF approximation in a "weak form" that allows us to study its 
stability by considering discrete versions of the W ,P -W ,q "duality pairing." This is equivalent 
to putting the QCF operator into a divergence form, which will indicate an interesting nonlocal 
effect of the atomistic to continuum interface. This nonlocal effect is the source of the lack of 
coercivity which we establish in Section HI based on the explicit construction of an unstable dis- 
placement. In Section [H we derive inf-sup stability results that are then combined, in Section [H 
with negative-norm truncation error estimates, to obtain optimal order error estimates for the 
force-based quasicontinuum approximation. We conclude, in Section [71 by showing the lack of a 
uniform inf-sup constant for all other common choices of duality pairings. 



We consider a one-dimensional atomistic chain whose 2M+1 atoms occupy the reference positions 
Xj = je, where e is the atomic spacing in the reference configuration, and which interact with their 
nearest and next-nearest neighbors. We denote the deformed positions by yj for j = — M, . . . , M. 
The boundary atoms are constrained by 



where F > is a macroscopic deformation gradient. The total energy of a deformation y 6 t 



2. The Force-Based Quasicontinuum Approximation 



y- M = -FMe and y M = FMe, 



is given by 
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where 



M 



E * 

j=-M+l 



Vi ~ Vj-i 



M 

+ E ^ 

j=-M+2 



Vi ~ V3-2 



(2.2) 



for a scaled two-body interatomic potential (ft (for example, the normalized Lennard-Jones potential 



(r) 



12^,-12 



e r 



2e r ) and external forces fj. The equilibrium equations are given by the force 



balance at the free atoms, 

i ? /(y) + /i=0 for j 

Vj = Fje for j 

where the atomistic force (per lattice spacing e) is given by 

ld£ a (y) if 



M + 1,...,M-1, 



(2.3) 



d Vi 



, / Vj+i - Vj 



+ 



I / Vj+2 - Vj 



i Vj- Vj-i 



+ 



i(Vj- Vj-2 



(2.4) 



In (|2.4p the undefined terms eft' \e~ l {y-M+i — V-M-i)) and cft'(e^ 1 (yM+i — Um-i)) are taken to be 
zero. 

We let Uj be a perturbation from the uniformly deformed state yj = Fje, that is, we define 

Uj = yj - Fje for j = -M, ... ,M. 

We linearize the atomistic equilibrium equations (|2.3|) about the deformed state y F , resulting in 

(L a u a )j = fj for j = -M + l,...,M -1, 



U : 



M, M, 



where (L a v)j, for a displacement v £ ' 
-Vj + i + 2vj - vj-i 



9 F 



for j : 
2M+1 , is given by 

V j+2 + Vj 



(2.5) 



+ 



>J 2F 



-Vj + i + 2vj - Vj-i 



-Vj + i + 2vj - Vj-i 



+ 4>2F 



-Vj+2 + 2Vj 



Vj-2 



+ 



'2F 



Vj - Vj-2 



-M + l, 

-M + 2,...,M-2, 
M - 1. 



Here and throughout we define <ft" F := 4>"{F) and (ft'^p '■= <ft"(2F), where (ft is the interatomic 
potential in (|2,2p . We assume that (ftp > and (ft'^p < 0, which holds for typical pair potentials 
such as the Lennard-Jones potential under physically relevant deformations. We remark that, for 



J 2F 



> 0, the system (|2.5p has a unique solution. This follows from (|3.ip and from Lemma 



(see also [9, 11] for an analysis of the periodic case which is similar). 

The local QC approximation uses the Cauchy-Born extrapolation rule to approximate the non- 
local atomistic model by a local continuum model [7,13,26,33]. In our context, this corresponds to 
approximating yj — yj-2 in (|2.2p by 2(yj — yj-i) and results in the local QC energy 



M 

E < 

j=-M+l 



Vj -Vj-i 



+ 



2 (yj - Vj-i) 



(2.6) 
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Note that the above expression has one more next-nearest neighbor term than f)2.2[) . This is because 
the atoms at j = —M + 1, M — 1 do not feel the effect of the boundary in the local approximation. 
The local quasicontinuum equilibrium equations are then given by 



if c (y) + /: 



3 

Vj 



for j = —M+ 1,...,M 

Fje for j = -M, M, 



1. 



where the local quasicontinuum force (per lattice spacing e) is given by 



if c (y) := 




, / Vj+i - Vj 



+ 2cl)' 



2(2/j+i ~ Vj) 



i yj- yj-i 



+ 2cf>' 



2 (Vj - Vj-i) 



(2.7) 



Linearizing the local quasicontinuum equilibrium equations (|2.7p about the deformed state y^ 
results in 

{L l q c u l qc) . = f . for i = _ M + l,...,M-l, 
for j = -M, M, 



uf = 



where (L lqc v)j, for a displacement v € M. 2M+1 , is given by 



{L lqc v) 



-v j+1 + 2 Vj - vj-x 



-M+ 1,...,M- 1. 



The increased efficiency of the local QC approximation is obtained when its equilibrium equa- 
tions (|2.7p are coarsened by reducing the degrees of freedom using piecewise linear interpolation 
between a subset of the atoms [7,26]. For the sake of simplicity of exposition, we do not treat 
coarsening in this paper. 

In order to combine the accuracy of the atomistic model with the efficiency of the local QC 
approximation, the force-based QC method decomposes the reference lattice into an atomistic 
region A and a continuum region C, and assigns forces to atoms according to the region they are 
located in. Since the local QC energy (j2.6H approximates i/j — Uj~2 in (12. ip by 2(yj — Vj-i), it is 
clear that the atomistic model should be retained wherever the strains are varying rapidly. The 
QCF operator is then given by [7, 8] 



Ff(y) 



and the QCF equilibrium equations by 



F 



gcf 



(y) + n 





Fje 



Ff{y) 



for j 
for j 



if j G A, 

if jeC, 



-M + l, 
-M, M. 



(2- 



,M-1, 



The force-based quasicontinuum approximation gets its name from the assignment of forces at 
the atoms in (|2.8p . Most other quasicontinuum approximations build a total energy by summing 
energy contributions from each region and compute forces on the atoms by differentiating the 
energy. However, F qc f is a non-conservative force field and cannot be derived from an energy [7]. 
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2.1. Artificial boundary conditions for the computational domain. For large atomistic 
systems, it is necessary to reduce the computational domain, even when using a coarse-graining 
method such as the QC approximation. The reduction of the computational domain requires the use 
of artificial boundary conditions to approximate the effect of the far field. The artificial boundary 
condition most commonly used in the QC method (and in computations using other atomistic to 
continuum approximations) sets the displacement to zero at the boundary of the computational 
domain, such as at the lateral boundary of the crystal in the nanoindentation problem reported 
in [18]. More accurate artificial boundary conditions such as given and analyzed in [21] do not seem 
to have yet been used in quasicontinuum computations. 

We chose to imitate the approach commonly used in the QC method, by choosing N <C M 
and < K < N — 1, letting {— N, . . . , N} be the computational domain, A = {— K, . . . , K} the 
atomistic region, and C = { — N + 1, . . . , N — 1} \ A the continuum region, and defining 



i7 c/ (y) + /i = for j = -N + l,...,N-l, 



(2.9) 



(2.10) 



Vj = Fje for j = -N, N, 

to be the QCF approximation on the computational domain. In this paper, we analyze the lin- 
earization of (|2.9I) . 

(L« c V c 0i = f ° r j = -N + l,...,N-l, 
uf f = u] for j = -N, N, 

where we have taken = N and u q ^ = ufj- so that we may ignore the error induced by the 
artificial boundary condition and exclusively focus on the error of the QC approximation. Note 
that, since atoms near the artificial boundary belong to C, only one boundary condition is required 
at each end. 

Setting e = 1/N throughout, we scale the problem [6] so that the size of the computational 
domain is of order 0(1). 

2.2. Notation. We use D : M. 2N+1 — ► M. 2N to denote the backward difference operator, defined by 

(Dv)j = Dvj = Vj ~ V J- X for j = -N + 1, . . . , N. 

We will frequently employ the weighted F-norms, 

/ N \ VP 

IMI^? '■= [ e Yl N P ' 1 - p < °°' 

V j=-N J 



max 

-N<j<N 



and the weighted inner product 



N 

(v,w) = ev i w r 

j=-N 

The definition of the difference operator D, of the norms ||v||^p and of the inner product (v,w) is 
extended, in an obvious way, for vectors v,w e M K , where K G N is arbitrary. For example, if v G 
R 2M+i then Dv G R 2M_ Moreover, in view of this convention, the higher order difference operators 
D 2 , etc., can be defined by successive application of D; for example, D 2 Vj = e 2 (vj — 2vj-i + Vj-2)- 
The subspace of M 2Ar+1 with homogeneous boundary conditions is denoted 

V = {v e R 2N+1 : V- N = v N = 0}. 



STABILITY, INSTABILITY, AND ERROR OF THE FORCE-BASED QC APPROXIMATION 



7 



For future reference we note that the following Poincare inequality holds [28, Lemma A. 3]: 

||v|^°o < |[|£>v[|^i for all v G V . (2.11) 

Furthermore, we note that the linear operator L qc f which has been defined above as a mapping 
from M 2JV+1 to M 2 ^ -1 will be considered below to be a mapping from R 2N+1 to V by the extension 

(L qcf v)- N = (L 9c/ v)jv = for v G R 2N+1 . 

With this in mind, (L qc fv,w) is well-defined for all v, w G IR 2Ar+1 . 

2.3. Pointwise consistency of the force-based QC approximation. The remarkable sim- 
plicity of the formulation of the force-based QC approximation is mirrored by its equally simple 
consistency analysis. Let u a be the solution to (j2.5f) (assuming <j>" F + 40 2 V > 0, this system is 
well-posed), then the truncation error t G M 2Ar+1 is defined by t-N = i/v = and 

tj = (L qcf u a - i)j = (L qcf u a - L a u a )j for j = -N + 1, . . . , N - 1, 

where L qc f\i a is understood by restricting u a to the computational domain. Since (L qc ^u a )j = 
(L a u a )j trivially holds for j G A we have tj = for j G A. For j G C, on the other hand, we have 



tj = (L qcf u a - L a u a )j = ^ F 



t -v% +1 + 2u« - -u a j+2 + 2u a j - u a j_ 2 



e (f> 2F 



e 2 e 2 



u a j+2 - ^u a j+l + Qu a j - Au°j_ x + u°j_ 2 



e 2 <P>i F {D A u a ) 3 , 



where (D 4 \r)j = (Z? 4 v)j + 2 is a fourth-order centered finite difference operator. Note also that u a 
is defined outside the computational domain. Thus, for p G [1, oo], we obtain an exact truncation 
error estimate, 

||t||, ? =e 2 |^ F |l|5 4 u a [| l P (c) , (2.12) 

where the label C indicates that the summation (or maximum) is only taken over atoms in the 
continuum region. 

We have presented (12.120 as a simple argument for the high accuracy of the QCF method, 
however, in the error analysis in Section [H] we will use a slightly sharper negative-norm estimate. 
We also note that it follows from the interior regularity of elliptic finite difference operators [34] 
that ||l) 4 u a ||^P( C ) is bounded in the continuum limit e — > if f is the restriction of a smooth function 
in a neighborhood of the continuum region C to the lattice points since the continuum region C is 
far from the boundary of the atomistic problem if N <g M. 

To estimate the error between the atomistic and QCF solution, we write 

L qc f(u a -u qc f)=t = 0(e 2 \^ F \). 

Hence, a uniform stability result for the operator L qc f in an appropriate norm would lead to an 
optimal error estimate. As we have already remarked in the introduction and will make precise 
in Theorem 14.11 L qc f is typically not coercive and we must therefore prove an inf-sup condition 
instead. To this end, we will factor the L qc f operator into divergence form, L qc f = D T E qc f D, 
where D is the discrete difference operator defined above. We will give conditions under which 
fiqcf j g row diagonally-dominant and which will lead to an inf-sup condition for L qc f . Interestingly, 
however, this approach only leads to uniform stability bounds if the t^°-t\ duality pairing is used, 
while the inf-sup constants for the £?-£f (1/p + 1/q = 1, l<p< oo) pairings are not uniform in N 
(cf. Section [7|). 
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Vw e V, 



for j 



0, 

-N, N, 



3. Divergence Form of the QCF Operator 
We will analyze the QCF equilibrium equations (|2.10p by putting them into a "weak form:" find 

u icf G JR2A+1 guch that 

(E qcf Du qcf , Dw) =(f , w) 
uf=u« 

where the linear operator E qcf : R 2N -> M 2Ar is chosen so that (E qcf Dv, L>w) = (L« c/ v, w) for all 
v G E 2iV+1 and w G V . We call £ gc/ the conjugate operator. This operator was previously derived 
for a Neumann problem in [7] and for a problem with mixed boundary conditions in [8]. 

To motivate the idea we briefly review the conjugate operator for the full atomistic model before 
deriving E qc * . The atomistic energy (and similarly all QC energies) can be written as functions of 
the strain Du, £ a (Du) := £ a (y), and its conjugate operator is defined by [7,8] 



ld£° 



e drj 



Thus, (E a (r))j is the negative of the force conjugate to the strain (Du)j. It follows from the chain 
rule that 



ld£ a (y) 
e 



-i 



d£ a 
dr j+1 

( li?(r)) 



d£ a 



Applying this calculation to the linearized operator L a , one can easily verify that (L a v,w) 
(E a Dv,Dw), for all v 6 R 2M+1 and w G R 2M+1 with W- M = % = 0, where 



E a 



y 2F 



1 



(3.1) 



From this representation we obtain immediately that, for v £ M 2A/+1 and for w € Vo, extended by 
zero outside the computational domain, the operator L a can be written in the "weak" form 

N 

(L a v, w) = WfDvj + (t>2 F {Dv j+1 + 2D Vj + Dvj-i)]Dwj. (3.2) 

This formula will be used in Section [6] to derive a negative-norm truncation error estimate. 

To find a representation for the QCF operator L qc f in terms of a conjugate operator we cannot 
simply carry out the same computation as above, even in the linearized case, since it is not related 
to any energy functional. Instead, we will first derive a "weak" form for L qc f from which it will be 
fairly straightforward to construct the conjugate operator. We begin by writing L qc f in the form 
L qc f = 4>" F L X + <% F L 2 , where 

,-2/ 



Vj+x + 2vj 



Uj-l 



)• 



-N + l,...,N -1, 



and 



4e 2 ( - Vj+i + 2vj - Vj-i) , 



e 2 ( - Vj+2 + 2vj 



0. 
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and deriving "weak" representations of the operators L\ and L 2 . 

Lemma 3.1. For all v € R 2Ar+1 and w G Vo the nearest neighbor and next-nearest neighbor 
interaction operators can be written in the form 

N 

(Liv,w)= eDvjDwj, and 

j=-N+l 

(L 2 v,w) = (L r 2 e9 v,w) + e 2 (D 3 v^ K+1 )w- K - e 2 (D 3 v K+2 )w K , 

where D 3 is the third-order backward finite difference operator, D 3 Vj = e~ 2 (Dvj — 2Dv 3 -\ -\-Dvj-2), 
and where L r 2 eg denotes the "regular" component of L2, 

-K K N 

{L r 2 e9 v,w) = ^2 e4DvjDwj+ ^ e(Dvj-i + 2Dvj + Dv j+1 )Dwj + ^ eADvjDwj. 

j=-N+l j=-K+l j=K+l 

Proof. We only prove the representation for L 2 . To simplify the notation, we will perform all 
manipulations only in the right half of the domain and indicate the remaining terms by dots, for 
example, 

K N-l 
(L 2 v, v) = ■ ■ ■ + ^ e(L 2 v) jWj + < L 2v)jWj. 

3=0 j=K+l 

The proof simply requires careful summation by parts, performed separately in the continuum 
and atomistic region. In the right half of the atomistic region, summation by parts yields 



3=0 



3=0 3=0 

K+2 K 

fVj - 

) W 3 

j=2 ' ' j=0 " 

K 

= . . . + e(D Vj + Dvj-i) (Dwj + Dwj-!) 
3=2 

- [(Dvk+i + Dv K )w K -i + {Dv K+2 + Dv K +i)w K ] 

K 

= h ^ e{Dvj+i + 2Dvj + Dvj-i)Dwj 

3=1 

- (Dv K +2 + 2Dv K+ i + Dv K )w K - 

Here, we also used the dots to indicate additional terms which would have canceled had we per- 
formed the calculation over the entire domain. A similar computation in the continuum region 
gives 

N-l N 

y~] €(L2'v)jVj = ^2 dDvjDwj + ADvk+iWk- 

j=K+l j=K+l 

Considering the symmetry of the problem, or by performing the same calculation in the left half of 
the domain, we obtain the stated result. □ 
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In order to find the conjugate operator, we only need to write wk and W-k in terms of the 
strains Dwj. This is achieved by connecting these displacements to the boundary, for example, we 
can use the identities 



N 



W K 



j=K+l 



and W-k 



-K 



eDwi 



-N+l 



Note, however, that there is no unique way of achieving this. Our choice above simply minimizes 
the number of non-zero entries for E qc f in each row, a fact that will become important later on. 
Thus, we obtain 



E qcf 



} 2F 



1 



1 



4 1-21 
5 -2 1 
1 2 1 
1 2 



With this choice we indeed obtain the identity (E qc f Dv, Dw) 
w G Vq. 



(3.3) 

1 

2 1 
-2 5 
-2 1 4 

-2 1 4 

(L 9 ^v, w) for all v G R 2N+1 and 



4. Lack of Coercivity 

Local minimizers of the atomistic energy are characterized, essentially, by the fact that the 
Hessian is positive definite. It can be shown that the coercivity of the operator L a in appropriate 
norms is independent of problem size N, provided that <j)" F + Afi^p > (see [9, 11] for a periodic 
problem; the proof for the Dirichlet boundary value problem is very similar). It can moreover 
be shown that the energy-based quasi-nonlocal approximation always inherits coercivity of the 
atomistic operator [11]. 

In Section El we will give conditions on (f) F and cft^p under which the force-based QC operator 
L qc f inherits stability (in a more general sense) uniformly in N, even though it does not inherit 
coercivity. In fact, as we show in the following theorem, the L qc f operator is not coercive whenever 
N is sufficiently large. That is, not only is it lacking uniform coercivity, but it is not even positive 
definite. This result clearly demonstrates why we need to work with a technically more involved 
inf-sup stability condition in our error analysis in later sections. 

Theorem 4.1. Suppose that <j)" F > and <\>'lp G R\ {0}; then, for sufficiently large N, the operator 
L qc f is not coercive. More precisely, there exist No G N, C\ > Ci > such that, for all N > Nq 
and 2 < K < N/2, 

-dN 1 ' 2 < inf (L 9c/ v,v) < -C 2 N 1 ' 2 . 
vsV 
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Proof. As in Section [3] we write L qc f = (j)" F L\ + (^'f-^ 2 - Since (L\v, v) = ||Z)v||^j, we need to 
concentrate on the next-nearest neighbor interaction operator L 2 . If we can show that L 2 is neither 
bounded above nor below, uniformly in N, and with the stated asymptotic behavior, then the 
upper bound follows. The lower bound follows from the fact that L\ is bounded while |(£ 2 v, v)| < 
CiN l / 2 \\Dv\\] 2 . Both of these facts are established in the following lemma. □ 

Lemma 4.1. Under the conditions of Theorem \4.1\ there exist positive constants c\, c 2 , independent 
of N , and lattice functions v + , v € Vo such that ||Dv + ||^2 = ||-Dv~||^ 2 = 1, 

(L 2 v+, v+) > ci(iV 1/2 - c 2 ), and (L 2 v~, v") < -ci(iV 1/2 - c 2 ). 

Moreover, these bounds are asymptotically optimal in that there exists a constant C3 > such that 

|(L 2 v,v)| < c 3 N 1/2 for allv € V with \\Dv\\g% = 1. 

Proof. We write (L 2 v, v) by setting w = v in Lemma [3. 11 The crucial observation is that the term 
vk(Dvk+2 — 2Dvk+i + Dvk) cannot be expressed as a quadratic form of strains supported at the 
interface, while all other terms are bounded in terms of (a constant multiple of) ||-Dv|| 2 2 . More 
precisely, we recall that 

(L 2 v, v) = (L r 2 eg v,v) - v K (Dv K+2 - 2Dv K+ i + Dv K ) + v- K (Dv- K -i ~ 2Dv_ K + Dv_ K+1 ), 
where K-L^Vi v )| — C 2||-D v || 2 2- Next, we construct the functions v 1 * 1 by choosing vk = 1 and so 
that the third difference in the bracket is of order N 1 ^ 2 . 

To this end, we set v = v + e l l 2 8K+i = v + A^ _1 / 2 <5a'+i, where 

f (N + j)/(N-K-2), j = - N ,...,-K-2 
Vj = { 1, j = -K-2,...,K + 2, 

( (N-j)/(N-K-2), j = K + 2,...,N, 

(that is, Vj = 1 in the atomistic region and the interface, and interpolates linearly between 1 and 
in the continuum region) and where $K+i,j = if j 7^ K + 1 and Sk+i,k+i = 1. In that case, 
1 1 Dv I \(2 is clearly uniformly bounded, and we obtain 

(L 2 v,v) = (L^ 9 v,v) +3iV 1 / 2 . 

Note that no terms at the left interface occur since v is a constant there. Upon appropriately 
rescaling by v + = v/[|JDv||^2 so that ||Dv + ||^2 = 1, we obtain 

(L 2 v+,v+) > -C2 + C1N 1 ' 2 . 

Setting v _ = c(v — e 1 ' 2 SK+i) gives the opposite bound. 

To prove the final statement, namely that these bounds are asymptotically sharp, we note that 
all terms of the type vxDvj are of order N 1 / 2 , 

\v K D Vj \ = e- 1/2 \v K \e 1/2 \Dvj\ < e~ 1/2 ||v||^o pv|^ < (2/e) 1/2 ||Z)v|||, 

where we used (I2.1ip and a weighted Cauchy-Schwartz inequality to bound [| v||^oc < \/2||Dv||_£2 . □ 

Remark 4.1. The proof of Lemma 14. II reveals that N needs to be of the order (I + I^/^fI) 2 before 
a loss of coercivity can occur. Although it may seem that this is typically a fairly large number, 
(1 + I^f/^'fI) 2 ^ s n °t so large for strains F near the edge of a stability region (such as near 
the critical strain at which the atomistic system "fractures" [5]), or more generally whenever the 
next-nearest neighbor interaction is not significantly dominated by the nearest neighbor interaction. 
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5. Stability of the Force-Based Quasicontinuum Solution 

We first recall a classical characterization of the norm of the inverse of an operator that we will 
use to prove the stability of the solution to the force-based quasicontinuum approximation. The 
proof is included for the sake of completeness. 

Lemma 5.1 (Inf-Sup Condition). Let W and V be finite dimensional normed linear spaces satis- 
fying dim W = dim V, and let L be a bounded linear operator from V to W where W is the dual 
of W. Suppose that 

inf sup (Lv,w) = 7 > 0. (5-1) 
„ vev weH / 

ll v llv-l ||w||w=l 

Then L is invertible and the solution u G V to Lu = f satisfies the stability bound 

1 

||u||y < — ||f||w where ||f||w := SU P (f 5 w )- 
7 wgVK 

ll w llw=i 

Proof. The inf-sup condition (|5.1|) implies that the nullspace of L must be trivial. Since a finite- 
dimensional linear operator between two spaces of the same dimension is invertible if and only if it 
is non-singular, we conclude that there is a unique solution u G V to Lu = f for every f G W'. 
If [|u||v = 0, then the stability bound is trivial. Otherwise, we have 

||f \\w' = SUp (Lu, w) = ||u||y SUp ( L ( .. .. ) ,W ) > 7||u||y. □ 
wgW wgW \ Vll u l|vV / 

||w||yy = l ||w||jy=l 

Next, we note that the range of the backward difference operator D is 



and therefore 



j=-N+l 



inf sup (L 9c/ v, w)= inf sup {E qcf Dv, Dw) 

vev we v vev we v 

||Dv|| £?0 =l ||£) W || 1= i \\Dv\\ e oo=l ||Dw||,i=l 



inf sup (E qc f£, rj). 

-1B>2JV _ m oiv 



H€llfS» =1 ||tj|Li=1 

The following lemma gives a bound on such an inf-sup constant, for a general matrix A. This 
result and its proof were inspired by [28, Sec. 3.1]. 

Lemma 5.2. Suppose that A G M_ 2Nx2N satisfies 

min (A u + A~A - max ^A± =: 7 > 0, 

where A~- = min(0, Ay) and Af- = max(0, Ay), then 

inf sup (A£, rj) > 7/2. 
ll£ll^ =1 [|t}||. 1 =1 
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Proof. Let £ G M^\{0} and choose p,q G {— JV+1, . . . , iV} such that £ p = maxj £j and £ g = min,- £j 



Since X^jL-Ar+i £j = 0' ^ follows that £ p > and £ g < 0. Moreover, let P = {j : £j > 0} and 
Q = {j : £■ < 0}. If we define 77 G by 



2e> 
_ J_ 
2e' 



0, 



i = q, 
otherwise, 



then 



j j 

>{4*&+E4& + E A ^}-{ A ^ + Y. A t^ + E 







jeP\{p} 






jeQ\{q} 






+ E 


{A«6 + E ^g}6 


+ E ^ 




j&Q 


jeP\{p] 




jeP 






A pp~ E I^pj 




161 + 


A? ~ E 


-Ei^-i]^ 




jeP\{p} 


jeP 




i6Q\{<?} 





> 7(1^1 + 1^1). 



□ 



From Lemma 15.21 and from (|3.3p , we can now deduce that 



inf sup 
veVo veV 
\\ Dv h°° =1 ||Dw|Li=i 



(L qcf v, w) > - 



min ((E^) u + ^(E« c %) - max^(i^)+ 



(5.2) 



Combining this estimate with Lemma 15 . 1 1 gives the following stability result. 



Theorem 5.1. Suppose that <j)" F + 80of > 0. Then the QCF system i2.10\) has a unique solution 
u qc f , which satisfies 



iDu^ll < 



2||f|| 



+ ■ 



+ 



N 



-N 



2N 



(5.3) 



where 



||f||* := sup (f,w). 

wGVq 

Proof. We write u qc f = u + u D where u G Vo and where = u a _ N + {u a N — u a _ N )(N + j)/(2N). 
Since u D is affine, it can be easily seen that L qc ^u D = 0. Hence, the system is equivalent to 
L qc f\i = f. In view of (|5.2p and Lemma 15.11 this has a unique solution, and we have the stability 
bound 



\Du 



qcf\ 



< ll-DulUoo + \\Du 



< 



2||f|L 



+ 



II 



N 



-N 



2N 



□ 
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6. Convergence 

The quasicontinuum error e qc f = u a — u 9 ^, where u a is again identified with its restriction to 
the computational domain whenever necessary, satisfies the equation 

L qc f e qcf ^j =tj, j = -N + 1, . . . , N - 1, 

(e^),- = 0, j = -N, N. 

Using (|2"T2"j) and (|2TL|) we see that the truncation error t = L qcf u a - f = (L qcf - L a )u a (but with 
tj\f = t_AT = 0) satisfies the negative norm estimate 

||t||* = SUp (t, W) < SUp |(t, w) = \ ||t||p = \<?\<\>2f\ 11-^ U<1 ||^(C)- 

weVo weVo 

||Dw|| f l=l jw|| £ oo=l 

However, we can get a slightly sharper result using the variational representations of the operators 
L a and L qc f derived in Section [3j 



Lemma 6.1. The truncation error satisfies the estimate 

||t||* <2e 2 |^ / F |p 3 u a || £?0(c ~ ) , 

where C = {-N + 2, . . . , -K + 1} U {K + 2, . . . , N + 1}. 

Proof. Using the "weak" forms of L a and L qc f derived in (13. 2h and in Lemma |3.1| we obtain 
<t,w) = <(L<^-L> a ,w) 

= <f>2 F l Yl e ( ~ + 2Du °i ~ Du j+l) Dw i + ( Du -k+i - 2Du a _ K + Du a _ K _ l )w^ K 

I j=-N+l 

N >| 

+ e (" + 2D«^ - + (- Du a K+2 + 2Du a K+1 - Du a K )w K \ 

j=K+l J 

< ^pVH^^dl^wll^ + 2||w|| <? >), 

where we used a weighted Holder inequality in the last step. Using (|2,lip to bound ||w||^oo we 
obtain the stated bound. □ 

Combining this negative-norm truncation error estimate with the stability estimate (|5.3p . we 
obtain the following result. 

Theorem 6.1. Suppose that <J) F + 84>2 F > 0. Then the atomistic problem \2. 5|) as well as the 
force-based quasicontinuum approximation $2.10\) have unique solutions, and they satisfy the error 
estimate 

\\D(u a - u qcf )\\ <4e 2 — — 'IhJtl 

4>F + ° ( t ) 2F 

As in Section [2] we note again that it follows from the interior regularity theory for elliptic finite 
difference operators [34] that IlL^u"!!^^ is bounded in the continuum limit e —* 0, provided that 

f is the restriction of a smooth function in a neighborhood of the continuum region C to the lattice 
points. 
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7. Estimates in Other Norms 



We conclude this paper by showing that our choice of norms with respect to which we analyzed 
the stability of the force-based QC approximation was, in some sense, unique. 

Theorem 7.1. Suppose that <f>" F > 0, <j)'l F £R \ {0}, and that 1 < p < oo, and 1 < q < oo so that 
| + | = 1. Then there exists a constant C > such that, for 2 < K < N/2, 

inf sup 

vev weVo 

\\Dv\\ e p=l ||d w ||^ = i 

Proof. We recall from Sections [3] and [5] that 



sup (L^v, w) = inf sup (E^t rj) < inf ||^ c/ ^||, P , 
l|Ov||, ? =i ||Dw|L 9 =i II€II#p=1||t,iu=i lieil#p=i 



inf 

vGVo 



where the second step follows from Holder's inequality. Therefore, we obtain the stated result from 
the following lemma. □ 



Lemma 7.1. Under the conditions of Theorem \ 7.1\ there exists a constant C > such that, for 
2<K< N/2, 



Proof. The terms causing this effect are the nonlocal terms extending from the atomistic to con- 
tinuum interface to the boundary. Hence, we define 

r 



to 



-1, 3 

-a, j 

0, j 
a, 3 

1, 3 



where a G K will be specified below, and £ 
for E qc f we see that 



-N + l,...,-K-l, 
-K, 

-K + l,...,K, 
K + l, 

K + 2,...,N, 

£/\\£\\e p - Recalling the matrix representation (|3.3 



E qcf i 





'-5 + 2a, 


3 


= -N + 1,...,-K -1 




-1 - 2a, 


3 


= -K, 




— a, 


3 


= -K + l, 


f€ + <t>2F < 


o, 


3 


= -K + 2,...,K -1, 




a, 


3 


= K, 




l + 2a, 


3 


= K + l, 




5 -2a, 


j 


= K + 2,...,N, 



from which we obtain 



|£ 9c/ £ll! = 2e(\a^ F \ p + \a<t>" F + (1 + 2a)^ F \ P 

+ (iV -ir-l)|^ + (5-2«y^|f). 
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A" +5(6" 

Choosing a = F 2 ^n 2F , and thereby canceling the term (N — K — 1)\4>" F + (5 — 2a)<p2 F \ p above, 
gives 

ll^ellj = 2e(|a0^r + + (1 + 2 «)^V| P ) • 

Moreover, since 

\\i\\ p £ p = 2e(N -K-1 + \a\ p ) > 2e(N/2 - 1 + |a| p ), 

we conclude that 

\p\ 1 /p 



• f ||^ c /,|| < / |«^r + |a^ + (l + 2a)^ F | P Y 
™2iv 11 *H<? - I iV/2-l + |a|f J 



fie 

ll€ll/?=i 



Conclusion 



We have presented a detailed stability and error analysis of the force-based QC method in one 
dimension. Although we were able to establish optimal order error estimates, we have also presented 
several "negative" results which are, in many respects, even more interesting. The present paper 
has focused exclusively on the force-based QC method, but we expect that the lack of coercivity 
(and more generally lack of stability in most norms) may be present in other force-based coupling 
methods such as [20,31] or the QM-MM coupling methods described in [4]. A careful study of these 
related methods is required to further understand and establish force-based coupling techniques as 
predictive tools in computational physics. 

Finally, let us remark on the the fact that we have only proven stability of the QCF method under 
the condition that qb" F + 8(j>2 F > 0. By contrast, the atomistic model is uniformly stable if, and only 
if 4>f + [Hi- Hence, we expect that our above condition is not sharp. Although we have 

established sharp characterizations of the stability of other QC methods in [11], we were unable 
to rigorously achieve the same for the force-based QC method as well. While our computational 
results reported in [11] do indicate that L qc > has only positive eigenvalues if <p" F + 4(/>2^ > 0, this 
is not enough to ensure stability, uniformly as N — > oo, in the regime <p" F + ^<p2 F > for the 
l^ 1 ' 00 -^ 1 ' 1 "duality pairing" or for any other norm with a recognizable continuum limit. 
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